/* Colluding Against Workers: Main model
	Delabastita & Rubens 
------------------------------------------*/

 
clear
set more off
cd "$station"

use ./data/temp/data_pf, clear 
set matsize 10000

xtset mineid yr 
xtdes

preserve
drop open
gen open = emp>0
drop if emp==.
collapse (sum) open, by(dq4id yr)
sum open, d
tab open
restore

 

sum uw* 
/* 1. Production function
--------------------------------*/

* First stage: agricultural wages and mining labor supply

sum seger* 
 
gen lsegers  = log(segers_nomwage_agr/cpi)
gen lsegers_ind  = log(segers_nomwage_ind/cpi)
  
preserve
collapse (mean)   lsegers* (sum) wl emp , by(yr)
tset yr
reg D.lsegers_ind D.lsegers


gen lemp = log(emp)
tset yr

gen Dlemp = D.lemp
gen Dlsegers = D.lsegers
gen Dlsegers_ind = D.lsegers_ind
 
reg D.lemp D.lsegers D.lsegers_ind  yr
 
reg Dlemp Dlsegers   , r

gen bsegers1 = _b[Dlsegers] 
gen sesegers1 = _se[Dlsegers] 
gen bsegersind1 = . 
gen sesegersind1 = .
local r2 = string(e(r2)) 		 
local r2_seg1 = substr("`r2'" ,1,4) 
local N_seg1 = e(N)

reg Dlemp Dlsegers Dlsegers_ind   , r

gen bsegers2 = _b[Dlsegers] 
gen sesegers2 = _se[Dlsegers] 
gen bsegersind2 = _b[Dlsegers_ind] 
gen sesegersind2 = _se[Dlsegers_ind] 
local r2 = string(e(r2)) 		 
local r2_seg2 = substr("`r2'" ,1,4) 
local N_seg2 = e(N)

gen seg = .
gen segind = . 
forvalues n = 1/2 {
replace seg = bsegers`n'
replace segind = bsegersind`n'
estpost su seg segind
est store bsegers`n'
replace seg = sesegers`n'
replace segind = sesegersind`n'
estpost su seg segind
est store sesegers`n'
}
 
label var seg "$\Delta$ log(Agricultural wage)"
label var segind "$\Delta$ log(Industrial wage)"

esttab bsegers1 sesegers1 bsegers2 sesegers2  using ./output/tab/table_pf_stage1.tex, replace  ///
mtitle("Est." "S.E." "Est." "S.E.") nolines prehead(  &  \multicolumn{2}{c}{ $\Delta$ log(Coal mining employment)  }  &  \multicolumn{2}{c}{ $\Delta$ log(Coal mining employment)  }    \\  )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs    ///
prefoot( && \\  R-squared  &	\multicolumn{2}{c}{`r2_seg1'} &	\multicolumn{2}{c}{`r2_seg2'}  \\  Observations &	\multicolumn{2}{c}{`N_seg1'}  &	\multicolumn{2}{c}{`N_seg2'} \\ &&\\  \hline   )  posthead( \hline  &&\\   ) 

 
reg Dlemp Dlsegers   , r
 
 
restore
   
  
** Specification 1: Cobb-Douglas   

xtset mineid yr
 
* OLS

reg lq lemp lwm   lk  
mat coef = e(b)
mat list coef
gen rho_ols1 = .
gen bl_ols1 = coef[1,1]
gen bm_ols1 = coef[1,2]
gen bk_ols1 = coef[1,3]
gen bc_ols1 = coef[1,4]
gen blk_ols1 = .
gen blm_ols1 = .

gen thl_ols1 = bl_ols1  
gen thm_ols1 = bm_ols1  
gen scale_ols1 = thl_ols1 +  thm_ols1 + bk_ols1

local r2 = string(e(r2)) 		 
local r2_ols1 = substr("`r2'" ,1,4) 
local N_ols1 = e(N)
 
* GMM
 
 gen lsegers_ag = log(segers_nomwage_agr/cpi)
 gen lsegers_ser = log(segers_nomwage_ser/cpi)
 
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm L.lsegers L.lk lk   )  conv_maxiter(50) //technique(bfgs)

 
 * sample sizes

di _N 	

gen dumobs = 0
replace dumobs = 1 if lq  ~=. & lemp ~=. & lwm~=.& lk~=.&lsegers ~=.
tab dumobs		

preserve 
keep if dumobs==1 
gen npf = 1
collapse (sum) npf, by(mineid)
local Nf_ols1 = _N
di `Nf_ols1'
restore 
 

replace dumobs = 0
replace dumobs = 1 if lq~=. & lemp ~=. & lwm~=.& lk~=.&lsegers ~=. &  L.lq ~=. & L.lemp ~=. & L.lwm~=.& L.lk~=.&L.lsegers ~=.
tab dumobs		
    
	
mat coef = e(b)
mat list coef
gen rho_bb1 = coef[1,1]
gen bl_bb1 = coef[1,2]
gen bm_bb1 = coef[1,3]
gen bk_bb1 = coef[1,4]
gen bc_bb1 = coef[1,5]
gen blm_bb1 = .
gen blk_bb1 = .

gen thl_bb1 = bl_bb1  
gen thm_bb1 = bm_bb1  

gen scale_bb1 = thl_bb1 +  thm_bb1 + bk_bb1 

* Number of firms   
   

preserve 
keep if dumobs==1 
gen npf = 1
collapse (sum) npf, by(mineid)
local Nf_gmm1 = _N
di `Nf_gmm1'
restore 
   
* R-squared
preserve 
xtset mineid yr
gen tfp = lq - bl_bb1*lemp-bk_bb1*lk - bc_bb1 -bm_bb1*lwm  
keep if tfp~=. & L.tfp~=.
gen qhat = tfp-lq
reg qhat lq 
egen avy = mean(lq)
egen tss = sum((lq-avy)^2)
egen rss = sum(tfp^2)
gen r2 = 1 - rss/tss
local r2 = string(r2) 		
di `r2'
local r2_gmm1 = substr("`r2'" ,1,4) 
restore
 

local N_gmm1 = e(N)
gen ltfp1 = lq - bl_bb1*lemp-bk_bb1*lk - bc_bb1 -bm_bb1*lwm 
gen tfp1 = exp(ltfp1)



* Model with depth control 

gen lsu =  s/u 
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)  -({bd})*(lsu - {rho}*L.lsu)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )      , ///
inst(  L.lemp  L.lwm  L.lsegers  lsu  L.lsu  L.lk lk  )  conv_maxiter(50) 

mat coef = e(b)
mat list coef
gen rho_bb12 = coef[1,1]
gen bl_bb12 = coef[1,2]
gen bm_bb12 = coef[1,3]
gen bk_bb12 = coef[1,4]
gen bd_bb12 = coef[1,5]
gen bc_bb12 = coef[1,6]

gen ltfp_bb12 = lq - bl_bb12*lemp-bk_bb12*lk - bc_bb12 -bm_bb12*lwm  - bd_bb12*lsu
gen tfp_bb12 = exp(ltfp_bb12)

gen nu_bb12 = ltfp_bb12 - rho_bb12 * L.ltfp_bb12
reg nu_bb12 L.nu_bb12 

* J statistic Overidentifying restrictions
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )      , ///
inst(  L.lemp  L.lwm  L.lsegers     L.lk lk  )  conv_maxiter(50) 
estat overid 

 
local j  = string(r(J)  ) 		 
local j_gmm1 = substr("`j'" ,1,4) 

local jp  = string(r(J_p)  ) 		 
local jp_gmm1 = substr("`jp'" ,1,4) 
 

di `r2'
  

** Fix returns to scale to nu = 1.05

gen nu = 1.05
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -(nu-{bl}-{bm})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk L.lsegers  )  conv_maxiter(50) 
   
local N_rts1 = e(N)
 local Nf_rts1 = `Nf_gmm1'

mat coef = e(b)
mat list coef
gen rho_rts1 = coef[1,1]
gen bl_rts1 = coef[1,2]
gen bm_rts1 = coef[1,3]
gen bk_rts1 = 1-bl_rts -bm_rts 
gen bc_rts1 = coef[1,4]
gen blm_rts1  = .
gen blk_rts1 = .

gen thl_rts1  = bl_rts  
gen thm_rts1  = bm_rts 
gen scale_rts1 = nu 

 estat overid 
 local j  = string(r(J)  ) 		 
local j_rts1 = substr("`j'" ,1,4) 

local jp  = string(r(J_p)  ) 		 
local jp_rts1 = substr("`jp'" ,1,4) 


 
 preserve 
xtset mineid yr
gen tfp = lq - bl_rts1*lemp-bk_rts1*lk - bc_rts1 -bm_rts1*lwm 
keep if tfp~=. & L.tfp~=.
gen qhat = tfp-lq
reg qhat lq 
egen avy = mean(lq)
egen tss = sum((lq-avy)^2)
egen rss = sum(tfp^2)
gen r2 = 1 - rss/tss
local r2 = string(r2) 		
di `r2'
local r2_rts1 = substr("`r2'" ,1,4) 
restore


   ** Store results to table

* cobb-douglas
 
foreach mod in "ols" "bb" "rts" {
foreach var in "bl" "bm" "bk"  "rho"   "thm" "thl" "scale"{
gen `var' = `var'_`mod'1
}
estpost su bl bm bk rho
est store pf_`mod'1
estpost su thl thm scale  
est store th_`mod'1
drop bl bm bk    thl thm scale rho
}

 

/* 2. Markups and markdowns
--------------------------------*/

gen am = wm/rev
gen al = wl/rev 
 
foreach mod in "bb" "ols" "rts" {
gen mu_`mod'1 = bm_`mod'1/am
gen md_`mod'1 = ((bl_`mod'1)/al)/mu_`mod'1
gen mp_`mod'1 = 1 + (mu_`mod'1-1)+(md_`mod'1-1) //md_bb1*mu_bb1
gen mdp_`mod'1 = (md_`mod'1 - 1)/md_`mod'1
gen mup_`mod'1 = (mu_`mod'1 - 1)/mu_`mod'1
}


 
sum mp_bb1, d

gen negmargin = mp_bb1 < 1 
replace negmargin = . if mp_bb1 == .
tab negmargin	// 12% 



sum mdp_bb1 


sum mp*, d 

foreach var of varlist am al {
egen `var'_p95 = pctile(`var'), p(95)
egen `var'_p05 = pctile(`var'), p(05)
}

sum mu_bb1, d 
sum mu_bb1 if am > am_p05 & am<am_p95 &   al > al_p05 & al<al_p95 , d
 
sum mu* md*, d
 
bys yr: egen agemp = sum(emp)
gen sempag = emp/agemp


foreach var in "mu" "md"  "mdp" "mup" {
foreach mod in "bb" "rts" "ols" {
forvalues m = 1/1 {
bys yr: egen  med`var'_`mod'`m' = median(`var'_`mod'`m')
bys yr: egen  av`var'_`mod'`m' = mean(`var'_`mod'`m')
 gen s`var'_`mod'`m' = `var'_`mod'`m'*sempag
bys yr: egen  ag`var'_`mod'`m' = sum(s`var'_`mod'`m')
}
}
}




 
* Wage comparison 
  
bys yr: egen agwl = sum(wl) 
gen aguwl = agwl/agemp
bys yr: egen avuwl = mean(wl/emp)
bys yr: egen avuwud = mean(uwud)

twoway connect scholliers yr || connect aguwl yr, yaxis(2) xline(1897, lwidth(medthick) lpattern(longdash) lcolor(orange_red))
  
preserve
collapse segers* *md* scholliers , by(yr)  

foreach var of varlist segers* *md* scholliers {
gen l`var' = log(`var')
}

bys yr: sum avmd_bb1
 
reg  lavmd_bb1 lsegers_nomwage_ind lsegers_nomwage_serv   yr
reg  lavmd_bb1 lsegers_nomwage_ind lsegers_nomwage_serv   yr
restore
 
sum mu* md*
 
 
/* Labor supply model */
 
gen degreve_imp_p  = degreve_imp_pq/degreve_imp_q
label var degreve_imp_p "Import price (BEF)"
twoway connect degreve_imp_p  yr, xline(1871 1875, lwidth(medthick) lpattern(longdash) lcolor(orange_red))	graphregion(color(white)) // demand shock due to international coal price hike in 1871-1875
*graph export ./output/fig/fig_demshock.pdf, replace
   
  
* Cournot model 

bys mineid: egen e_union = max(dunion)
bys mineid: egen e_car = max(dcar)
  
forvalues m = 1/4 {
preserve
collapse(sum) emp wl drail dtram e_union dcar e_car (mean) cpi segers*  degreve*  , by(dq`m'id   yr)
gen p98 = yr>=1898
gen dcar_p98 = dcar*p98
gen e_car_p98 = e_car*p98
gen demshock = yr>= 1871 & yr<= 1875
gen lemp = log(emp)
gen luw = log(wl/emp)
gen lqimp = log(degreve_imp_q)
gen pimp = (degreve_imp_p/cpi)
gen lpimp = log(pimp)
gen lsegers  = log(segers_nomwage_agr/cpi)
gen lcpi = log(cpi)  
 
reg lqimp lpimp
predict demshock2, resid

gen ptrans =  ((degreve_trans_pq/degreve_trans_q)/cpi)
gen lptrans = log(ptrans)

reg luw lemp e_car p98 demshock e_car_p98
local becar_ols_dq`m' = _b[e_car]
local bp98_ols_dq`m' = _b[p98]
local bw_ols_dq`m' = _b[lemp]
local sew_ols_dq`m' = _se[lemp]

* sample sizes

di _N 	
gen dumobs = 0
replace dumobs = 1 if lemp  ~=. & luw ~=. & e_car ~=. & p98 ~=. & demshock  ~=. & e_car_p98 ~=.  
tab dumobs		
 

ivregress 2sls luw (lemp =  demshock  e_car_p98  )  e_car p98    

* driscoll kraay standard errors

xtset dq`m'id yr 
ivreg2 luw (lemp = demshock e_car_p98) e_car p98  , dkraay(2)  r

local N_ls_dq`m' = e(N)
 
local becar_iv_dq`m' = _b[e_car]
local bp98_iv_dq`m' = _b[p98]
local bw_iv_dq`m' = _b[lemp]
local sew_iv_dq`m' = _se[lemp]

mat jmat = e(j)
mat jmat2 = e(jp)

mat list jmat

local j = string(jmat[1,1])
local jp  = string(jmat2[1,1]) 		 
 
local j_iv_dq`m' = substr("`j'" ,1,4) 
local jp_iv_dq`m' = substr("`jp'" ,1,4) 
 

 
reg lemp demshock e_car_p98 e_car p98, r
local F  = string(e(F))
local F_ls_dq`m' = substr("`F'",1,3)





 
restore
gen bw_ols_dq`m' = `bw_ols_dq`m''
gen sew_ols_dq`m' = `sew_ols_dq`m''

gen bw_iv_dq`m' = `bw_iv_dq`m''
gen sew_iv_dq`m' = `sew_iv_dq`m''


}

 di `j_iv_dq4'
 di `jp_iv_dq4' 
 
* markdowns
 
forvalues m = 1/4 {
bys dq`m'id yr: egen emp_dq`m' = sum(emp)
gen semp_dq`m' = emp/emp_dq`m'
gen psic_dq`m' = 1+bw_iv_dq`m'*semp_dq`m'
gen psim_dq`m' = 1+bw_iv_dq`m'
gen kappa_dq`m' = (md_bb1- psic_dq`m'    )/( psim_dq`m' - psic_dq`m') 
gen kappa_rts_dq`m' = (md_rts1- psic_dq`m'    )/( psim_dq`m' - psic_dq`m') 
foreach var in "psic" "psim" "kappa" "kappa_rts" {
bys yr: egen av`var'_dq`m' = mean(`var'_dq`m')
bys yr: egen med`var'_dq`m' = median(`var'_dq`m')
}
}
 
 
* market- ad firm-level labor supply elasticities:
 
gen supelast_ols_m = 1/(bw_ols_dq4)
gen supelast_ols_f = 1/(bw_ols_dq4*semp_dq4)

gen supelast_iv_m = 1/(bw_iv_dq4)
gen supelast_iv_f = 1/(bw_iv_dq4*semp_dq4)

egen avsupelast_iv_f = mean(supelast_iv_f)
local avupelast_iv_f = avsupelast_iv_f

egen avsupelast_ols_f = mean(supelast_ols_f)
local avupelast_ols_f = avsupelast_ols_f
 
di `F_ls_dq4'
gen bw = bw_ols_dq4 
estpost su bw 
est store labsup_ols
replace bw = sew_ols_dq4 
estpost su bw
est store labsup_ols_se

replace bw = bw_iv_dq4 
estpost su bw 
est store labsup_iv
replace bw = sew_iv_dq4 
estpost su bw
est store labsup_iv_se



* locals for labor supply table
local supelast_ols_f = string(avsupelast_ols_f)
local supelast_ols_f_short = substr("`supelast_ols_f'",1,6)
local supelast_iv_f = string(avsupelast_iv_f)
local supelast_iv_f_short = substr("`supelast_iv_f'",1,6)
 
  

sum medpsi* medkappa*
sum bl*

sort yr 


 
  
* Evolution of markdowns 
  
gen lkappa_dq4 = log(kappa_dq4+1) 

sum lkappa* 
 
 
forvalues m = 1/4 { 

label var medpsic_dq`m' "Non-collusive markdown"
label var medpsim_dq`m' "Fully-collusive markdown"
label var medmd_bb1 "Actual markdown"
label var avpsic_dq`m' "Non-collusive markdown"
label var avpsim_dq`m' "Fully-collusive markdown"
label var avmd_bb1 "Actual markdown"

twoway connect medpsic_dq`m' medpsim_dq`m' medmd_bb1 yr, lpattern(solid solid solid) mlcolor(navy maroon green) mfcolor(white white white) msize(medium msmall medium) msymbol(circle diamond square)  graphregion(color(white))  xline(1897, lpattern(longdash)) ///
 xscale(titlegap(*5)) ylabel(,angle(horizontal)) ytitle("Median markdown") xtitle("Year")
 graph export ./output/fig/fig_mmd_decomp_dq`m'.pdf, replace
 
 twoway connect avpsic_dq`m' avpsim_dq`m' avmd_bb1 yr, lpattern(solid solid solid) mlcolor(navy maroon green) mfcolor(white white white) msize(medium msmall medium) msymbol(circle diamond square)  graphregion(color(white)) ///
 xscale(titlegap(*5)) ylabel(,angle(horizontal))
* graph export ./oligopsony_ir/fig/fig_amd_decomp_dq`m'.pdf, replace
 
}



  

  
* Store markup moments

 

foreach mod in "bb" "ols" "rts" {
gen medmd = medmd_`mod'1 
gen avmd = avmd_`mod'1
gen medmu = medmu_`mod'1 
gen avmu = avmu_`mod'1 

estpost su medmd avmd medmu avmu 
est store md1_`mod'
drop medmd avmd medmu avmu 
}

 
 
 
 
/* D. Drivers of monopsony power */

** D1. Market concentration and integration 

preserve 
 drop open
gen open= q>0
replace open = . if q==.
drop if open==. | open==0
 
forvalues m = 1(1) 4 {
bys dq`m'id yr: egen n_dq`m' = sum(open)
tab n_dq`m' , gen(nm_dq`m'_)
}

gen urb = dq2=="40" | dq2=="50" 
bys mineid: egen urban = max(urb)

foreach var of varlist md* mu* {
gen l`var' = log(`var')
}
 
reg lmd_bb1 nm_dq3_1 nm_dq3_2 nm_dq3_3 drail dtram urban  i.yr, vce(cluster dq3id )
gen th_rail1 = _b[drail]
gen th_tram1 = _b[dtram]
gen th_nm11 = _b[nm_dq3_1]
gen th_nm21 = _b[nm_dq3_2]
gen th_nm31 = _b[nm_dq3_3]
gen th_urban1 = _b[urban]
local N_conc1 = e(N)
local r2  = string(e(r2))
local r2_conc1 = substr("`r2'",1,4)

xtreg lmd_bb1 nm_dq3_1 nm_dq3_2 nm_dq3_3 drail dtram urban   i.yr, fe r


gen th_rail2 = _b[drail]
gen th_tram2 = _b[dtram]
gen th_nm12 = _b[nm_dq3_1]
gen th_nm22 = _b[nm_dq3_2]
gen th_nm32 = _b[nm_dq3_3]
gen th_urban2 = _b[urban]
local N_conc2 = e(N)
local r2  = string(e(r2))
local r2_conc2 = substr("`r2'",1,4)
 
* Store estimates
forvalues m = 1/2 {
rename (th_rail`m'  th_tram`m'  th_urban`m'  th_nm1`m'  th_nm2`m'  th_nm3`m' ) ( th_rail  th_tram  th_urban  th_nm1 th_nm2 th_nm3) 
estpost su  th_rail   th_tram  th_urban  th_nm1 th_nm2 th_nm3 
est store conc`m' 
drop th_rail  th_tram  th_urban  th_nm1 th_nm2 th_nm3 
}
restore 
 
** Collusion: cross-sectional variation

 
preserve
 tab e_union 
 bys e_union : egen qunion = sum(q)	
 egen qa = sum(q)
 gen squnion = qunion/qa
 bys e_union: sum squnion 
restore 
 
 drop e_union
   bys mineid: egen e_union=max(dunion)
 
 
gen lmd_bb1 = log(md_bb1)
reg lmd_bb1 e_union dcar yr, r

gen th_union1 = _b[e_union]
gen th_car1 = _b[dcar]
local N_col1 = e(N)
local r2  = string(e(r2))
local r2_col1 = substr("`r2'",1,4)

reg lmd_bb1  e_union dcar i.yr , r

gen th_union2 = _b[e_union]
gen th_car2 = _b[dcar]
local N_col2 = e(N)
local r2  = string(e(r2))
local r2_col2 = substr("`r2'",1,4)

bys dunion yr: egen medmd_bb1u = median(md_bb1)
bys dunion yr: egen avmd_bb1u = mean(md_bb1)

twoway connect medmd_bb1u yr if dunion==0 || connect medmd_bb1u yr if dunion==1, yline(1)
twoway connect avmd_bb1u yr if dunion==0 || connect avmd_bb1u yr if dunion==1, yline(1)
  
  
** Markdown-size correlation

gen lsemp_dq4 = log(semp_dq4)
egen mid = group(dq4id yr)

* (a) all firms

reg lmd_bb1 lsemp_dq4  yr
gen mdcorr1a = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr1a = substr("`r2'" ,1,4) 
local n_mdcorr1a = e(N)

areg lmd_bb1 lsemp_dq4 yr, absorb(dq4id)
gen mdcorr1b = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr1b = substr("`r2'" ,1,4) 
local n_mdcorr1b = e(N)

areg lmd_bb1 lsemp_dq4 yr, absorb(mid)
gen mdcorr1c = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr1c = substr("`r2'" ,1,4) 
local n_mdcorr1c = e(N)

* (b) non-cartel firms

reg lmd_bb1 lsemp_dq4  yr if dcar==0
gen mdcorr2a = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr2a = substr("`r2'" ,1,4) 
local n_mdcorr2a = e(N)

areg lmd_bb1 lsemp_dq4 yr if dcar==0, absorb(dq4id)
gen mdcorr2b = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr2b = substr("`r2'" ,1,4) 
local n_mdcorr2b = e(N)

areg lmd_bb1 lsemp_dq4 yr if dcar==0, absorb(mid)
gen mdcorr2c = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr2c = substr("`r2'" ,1,4) 
local n_mdcorr2c = e(N)

* (c) cartel firms

reg lmd_bb1 lsemp_dq4  yr if dcar==1
gen mdcorr3a = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr3a = substr("`r2'" ,1,4) 
local n_mdcorr3a = e(N)

areg lmd_bb1 lsemp_dq4 yr if dcar==1, absorb(dq4id)
gen mdcorr3b = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr3b = substr("`r2'" ,1,4) 
local n_mdcorr3b = e(N)

areg lmd_bb1 lsemp_dq4 yr if dcar==1, absorb(mid)
gen mdcorr3c = _b[lsemp_dq4]
local r2 = string(e(r2)) 		 
local r2_mdcorr3c = substr("`r2'" ,1,4) 
local n_mdcorr3c = e(N)

 

foreach n in "a" "b" "c" {
forvalues m = 1/3 {
gen mdcorr = mdcorr`m'`n'
estpost su mdcorr 
est store mdcorr_`m'`n'
drop mdcorr 
}
} 



 
** Collusion: changes over time

replace dumobs = 0
replace dumobs = 1 if lmd_bb1 ~=. 
tab dumobs	

replace dumobs = 0
replace dumobs = 1 if lmd_bb1 ~=. & e_union ~=. & dcar ~=.
tab dumobs		

 
*reg lmd_bb1 yr ,  r
reg lmd_bb1  yr ,  r
local N_time1 = e(N)
local r2  = string(e(r2))
local r2_time1 = substr("`r2'",1,4)
gen th_lyr = _b[yr] 


estpost su th_lyr
est store time1


gen dyr1 = yr<1865 & yr>=1855
gen dyr2 = yr<1875 & yr>=1865
gen dyr3 = yr<1885 & yr>=1875
gen dyr4 = yr<1895 & yr>=1885
gen dyr5 = yr<1905 & yr>=1895
gen dyr6 = yr<1915 & yr>=1905

reg lmd_bb1  dyr1 dyr2 dyr3 dyr4 dyr5 dyr6 
local N_time2 = e(N)
local r2  = string(e(r2))
local r2_time2 = substr("`r2'",1,4)
gen th_yr1 = _b[dyr1] 
gen th_yr2 = _b[dyr2] 
gen th_yr3 = _b[dyr3] 
gen th_yr4 = _b[dyr4] 
gen th_yr5 = _b[dyr5] 
gen th_yr6 = _b[dyr6] 

estpost su th_yr*
est store time2


** Collusion and employer associations: comparison pre- and post-cartel period:

reg lmd_bb1  e_union  i.yr i.dq1id if yr<1898 & dcar~=., r

gen th_union3 = _b[e_union]
local N_col3 = e(N)
local r2  = string(e(r2))
local r2_col3 = substr("`r2'",1,4)

reg lmd_bb1   e_union  i.yr i.dq1id if yr>=1898 & dcar~=., r
gen th_union4 = _b[e_union]
local N_col4 = e(N)
local r2  = string(e(r2))
local r2_col4 = substr("`r2'",1,4)

* Store estimates

forvalues m = 1/2 {
rename (th_union`m'  th_car`m'   ) ( th_union th_car) 
estpost su  th_union th_car
est store col`m' 
drop th_union th_car
}

forvalues m = 3/4 {
rename (th_union`m'      ) ( th_union  ) 
estpost su  th_union  
est store col`m' 
drop th_union  
}



* D4. Factor bias?
 gen lkemp = log(k/emp)
 gen lhp = log(hp  )
 
 corr md_bb1 k_ventilation_hp	  
 corr md_bb1 k_extraction_hp	   
 corr md_bb1 k_excavation_hp	  
 
   
  xtreg lk  lmd_bb1  i.yr, fe r
xtreg lhp   lmd_bb1 i.yr, fe r


* D3. Politics
 
gen dsoc  = ssoc1899 >=0.50
replace dsoc  = . if ssoc1899 == . 
 
bys dsoc  yr: egen avkappa1_dsoc = mean(kappa_dq4)
bys dsoc  yr: egen medkappa1_dsoc = median(kappa_dq4)
 
 twoway connect medkappa1_dsoc yr if dsoc==0, lwidth(medium) msize(medium) mcolor(white) mlcolor(navy) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) lcolor(navy) || ///
 	connect medkappa1_dsoc yr if dsoc==1, lwidth(medium) msize(medsmall) msymbol(diamond) mcolor(white) mlcolor(maroon) ///
 	xline(1894 1899, lwidth(medthick) lpattern(longdash) lcolor(orange_red)) legend(label(1 "Other majority") label(2 "Socialist majority")) ///
 	graphregion(color(white)) legend(col(2)) xscale(titlegap(*5)) ytitle("Collusion") xtitle("Year")
 graph export ./output/fig/fig_pol.pdf, replace

 twoway connect avkappa1_dsoc yr if dsoc==0, lwidth(medium) msize(medium) mcolor(white) mlcolor(navy) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) lcolor(navy) || ///
 	connect avkappa1_dsoc yr if dsoc==1, lwidth(medium) msize(medsmall) msymbol(diamond) mcolor(white) mlcolor(maroon) ///
 	xline(1894 1899, lwidth(medthick) lpattern(longdash) lcolor(orange_red)) legend(label(1 "Other majority") label(2 "Socialist majority")) ///
 	graphregion(color(white)) legend(col(2)) xscale(titlegap(*5)) ytitle("Collusion") xtitle("Year")


 
* E. Variance decomposition

sum dq4id

gen uw = wl/emp 

gen luw = log(wl/emp)

gen luwud = log(uwud)

reg luwud i.yr 
gen r2_w1 = e(r2)
reg luwud i.yr i.dq4id  
gen r2_w2 = e(r2)
areg luwud i.yr i.dq4id  , absorb(mineid)
gen r2_w3 = e(r2)

areg luwud i.yr i.dq4id  , absorb(mineid)

 
xtset mineid yr

*drop sumvar avvar bfe sd* var* sumsd* sdw varw pctvar*

reg luwud i.mineid    , r

mat b = e(b)
mat define btrans = b'
mat V = e(V)
mat list V 
mat list b

* Average standard error of the firm fixed effects
gen sumvar = 0 
local nfe = `= rowsof(V)'-1
di `nfe'	// number of fixed effects + 1 

forvalues n = 1/`nfe' {

replace sumvar = sumvar + (V[`n',`n'])
}
gen avvar = sumvar/`nfe'

* Variance of the firm fixed effects 

svmat btrans, names(bfe)
replace bfe = . if _n>`nfe'
replace bfe = . if _n==1
brow bfe 
egen sdfe = sd(bfe) 
gen varfe = sdfe^2	// "naive" variance of firm fixed effects
sum varfe
 
gen sdfe_bc = sqrt(varfe - avvar)
sum sdfe_bc

gen sumsd = sqrt(sumvar)
gen sumsd_cor = sqrt(sumvar-avvar)
sum sumsd sumsd_cor 
sum sumvar 
 
egen sdw = sd(luwud)
gen varw = sdw^2
sum varw  

gen pctvarw = (varfe - avvar )/(varw ) 
sum pctvarw 
 
 

reg lmd_bb1 i.yr
gen r2_m1 = e(r2)
reg lmd_bb1 i.yr i.dq4id  
gen r2_m2 = e(r2)
areg lmd_bb1 i.yr   , absorb(mineid)
gen r2_m3 = e(r2)


forvalues m = 1/3 {
gen r2_w  = r2_w`m' 
gen r2_m  = r2_m`m' 
estpost su r2_w r2_m
est store r2_`m'
drop r2_w r2_m
}
   
 
   
* Observed exit

xtset mineid yr

gen dx =  emp ~=. & F.emp==. 
replace dx = . if emp==.
replace dx = . if yr==1913	// last year

replace dx =  emp ~=. & F2.emp==.  if yr==1849  
replace dx =  emp ~=. & F5.emp==.  if yr==1885   
replace dx =  emp ~=. & F3.emp==.  if yr==1899   

sum dx 
bys yr: sum dx 
bys yr: egen avdx  = mean(dx )

* weight exit by year
bys yr: egen empyr = sum(emp)
gen sempyr = emp/empyr

gen wdx = dx*sempyr
bys yr: egen agdx  = sum(wdx)
bys yr: egen check = sum(sempyr)

label var avdx "Annual exit share of firms"
label var agdx "Annual exit share of employment"

twoway connect avdx agdx yr, xline(1898, lwidth(medthick) lpattern(longdash) lcolor(orange_red)) graphregion(color(white)) msymbol(circle diamond) msize(medium msmall) mcolor(white white) mlcolor(navy maroon) lcolor(navy maroon) ylabel(,angle(horizontal)) xtitle("Year") ytitle("Exit rate (Annual)") xscale(titlegap(*5))
graph export ./output/fig/fig_exityr.pdf, replace  

* observed exit smoothed over 4 year time horizon 


gen t = round(yr/4)*4
brow yr t   

bys t : egen avdx_t = mean(dx)
bys t : egen agdx_t = mean(agdx)

label var avdx_t "Exit share of firms"
label var agdx_t "Exit share of employment"

twoway connect avdx_t agdx_t t, xline(1898, lwidth(medthick) lpattern(longdash) lcolor(orange_red))  graphregion(color(white)) msymbol(circle diamond) msize(medium msmall) mcolor(white white) mlcolor(navy maroon) lcolor(navy maroon) ylabel(,angle(horizontal)) xtitle("Year") ytitle("Exit rate (Four-year window)") xscale(titlegap(*5))
graph export ./output/fig/fig_exitt.pdf, replace  

	
	
* Observed variable profits 

gen vp  = rev - wl-wm 
sum vp , d 
gen negvp  = vp <0
replace negvp  = . if vp ==.
tab negvp 

 
** Counterfactual exercise
 
/* how much does symmetric firms assumption matter? */

 

drop nf 
gen active = q>0 
bys dq4id yr: egen nf = sum(active)
gen s2 = semp_dq4^2
bys dq4id yr: egen hhi = sum(s2)
 
gen agmd_sym = (psim_dq4-1)/nf + 1  
gen agmd_asym = (psim_dq4-1)*hhi + 1  

preserve 
collapse(mean) agmd_sym agmd_asym, by(dq4id yr)
sum agmd_sym agmd_asym* 

twoway kdensity agmd_sym ||kdensity agmd_asym

sum agmd_sym agmd_asym, d
restore 
 
  
/* CF1: Exogenous P: */


save ./data/temp/data_cf_temp, replace


*preserve    
use ./data/temp/data_cf_temp, clear
  
drop if open==0|open==.
 
 gen bl = bl_bb1 
gen bm = bm_bb1
 gen psi = bw_iv_dq4

 collapse(sum) open rev wl   emp   (mean) winv bl bm psi (median) kappa_dq4 , by(dq4id yr)  
 
 rename open n
  
 gen uw = wl/emp
 gen uwbar = uw /(emp)^psi 
 
gen theta0 = 1/n	// cournot
gen theta1 = 1		// cartel
gen theta2 = 1/n + (1-1/n)*0.4	// pre-cartel collusion
gen theta3 = 1/n + (1-1/n)*kappa_dq4 	// estimated theta 
gen theta4 = 1/(n+1) + (1-1/(n+1))*kappa_dq4  // theta if there would be one additional entrant in the market

gen wm = bm* rev 

forvalues m = 0/4  {
gen uw`m' = (uwbar)^(1/(1+psi)) *((bl*rev)/(1+psi*theta`m'))^(psi/(1+psi))
gen emp`m' = (uw`m'/uwbar)^(1/psi)
gen vp`m' = rev-uw`m'*emp`m' - wm  
}


* estimate fixed costs 
 
gen fc_upper = vp3/n  
gen fc_below = vp4/n

gen fc_mid = (fc_upper + fc_below)/2
sum fc_mid, d 

sum winv, d 
sum winv if fc_mid~=., d

gen fcm_upper = fc_upper/1000000
gen fcm_below = fc_below/1000000

egen medfcm_upper = median(fcm_upper)
egen medfcm_below = median(fcm_below)

estpost su medfcm_upper medfcm_below 
est store fc_ex

gen lfc_mid = log(fc_mid)
gen lwinv = log(winv)

sum fc_mid winv 
reg lfc_mid lwinv
corr lfc_mid lwinv

  

label var lfc_mid "log(Estimated fixed cost)" 
label var lwinv "log(Capital investment)" 
gen x45 = lfc_mid 
label var x45 "45 degree line" 

twoway scatter lfc_mid lwinv||line lfc_mid x45, graphregion(color(white)) mfcolor(white) xtitle("log(Capital investment)") legend(off) xscale(titlegap(*5)) ylabel(,angle(horizontal)) 
graph export ./output/fig/fig_fc_winv_exg.pdf, replace  
 

* exit probabilities
forvalues n = 0/4 {
gen dx`n' = fc_mid*n > vp`n'
replace dx`n' = . if fc_mid==.|vp`n'==.
}
sum dx* 					 

 
 
* save exit probabilities to table
 
foreach var of varlist dx* {
egen av`var'=mean(`var')	
}

 
gen ratx01 = (avdx0/avdx1  - 1) 
  
gen ratx21 = (avdx2/avdx1 - 1) 
 
foreach var in "ratx"  {
gen `var' = `var'01 	
}

 

estpost su ratx 
est store ratx01_ex 

foreach var in "ratx"  {
replace `var' = `var'21	
}

estpost su ratx 
est store ratx21_ex 



 
 
foreach var of varlist uw0 uw1 uw2 {
bys yr: egen av`var' = mean(`var')
}

twoway connect avuw0 avuw1 avuw2  yr, xline(1898)

sum uw*
sum emp*

gen ratw10 = uw1/uw0 -1	// collusion vs. cournot
gen ratl10 = emp1/emp0 -1	//  

gen ratw12 = uw1/uw2 -1	// collusion vs. pre-cartel
gen ratl12 = emp1/emp2 -1	//  

sum rat*

replace uw1 = uw2 if yr<1898
drop avuw* 
  
  foreach var of varlist uw1 uw2 {
bys yr: egen av`var' = mean(`var')
}
  
* add counterfactuals to table

foreach var in "ratw" "ratl"  {
gen `var' = `var'10	
}

estpost su ratw ratl  
est store cf10_ex 

foreach var in "ratw" "ratl"  {
replace `var' = `var'12	
}

estpost su ratw ratl  
est store cf12_ex 

label var ratw "Relative wage change"
label var ratl "Relative employment change"
  
esttab cf10_ex cf12_ex  using ./output/tab/table_cf.tex, replace  ///
mtitle("Cournot" "Pre-1898 conduct"  ) nolines prehead(\textit{Panel A: Exogenous price }    &  \multicolumn{2}{c}{Comparison of cartel to: }      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&   ///
\\     \hline    )  posthead( \hline  &&&\\   ) 
 
twoway connect avuw2 avuw1   yr, xline(1898)

collapse(mean) ratw10 ratw12 ratl10 ratl12 fc*, by(yr dq4id)
save ./data/temp/ratw_exg, replace 
 
*** Endogenous P

* Counterfactuals  
 
use ./data/temp/data_cf_temp, clear
  
drop if open==0|open==.
 
 gen bl = bl_bb1 
 gen bm = bm_bb1 
 gen bk = bk_bb1 
gen psi = bw_iv_dq4

  
  gen p98 = yr>=1898

  gen e_car_p98 = e_car*p98
 

 
 
 
collapse(sum)  k_ventilation_hp open rev q wl emp (mean) tfp1 bl bm bk psi winv (median)   kappa_dq4 , by(dq4id yr)  


gen uw = wl/emp
gen p = rev/q

gen lp =log(p)
gen lq = log(q)

gen lk_ven = log(k_ventilation_hp)

reg lp lq 
gen eta_ols = _b[lq]
gen eta_ols_se = _se[lq]

gen ltfp1 = log(tfp1)

reg lq ltfp1 , r
local F  = string(e(F))
local F_coaldemand = substr("`F'",1,4)
local N_coaldemand_ols = e(N)

 
ivregress 2sls lp (lq = ltfp1) yr, r


gen eta_iv = _b[lq]
gen eta_iv_se = _se[lq]

local N_coaldemand_iv = e(N)

gen eta = eta_ols 
estpost su eta 
est store eta_ols 

replace eta = eta_ols_se 
estpost su eta 
est store eta_ols_se

replace eta = eta_iv 
estpost su eta 
est store eta_iv

replace eta = eta_iv_se
estpost su eta 
est store eta_iv_se 

replace eta = eta_iv


* add results to table 




 label var eta "log(Output)"
esttab eta_ols eta_ols_se eta_iv eta_iv_se  using ./output/tab/table_coaldemand.tex, replace ///
mtitle("Est." "S.E." "Est." "S.E.") nolines prehead(   &  \multicolumn{2}{c}{ log(Price)}   &  \multicolumn{2}{c}{ log(Price)}      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&\\ Method &	\multicolumn{2}{c}{OLS} &	\multicolumn{2}{c}{IV}  ///
\\  First-stage F-statistic && &	\multicolumn{2}{c}{`F_coaldemand'}\\   Observations &	\multicolumn{2}{c}{`N_coaldemand_ols'} &	\multicolumn{2}{c}{`N_coaldemand_iv'}  ///
\\   \hline    )  posthead( \hline  &&&\\   ) 
 



rename open n


gen uwbar = uw /(emp)^psi 
gen pbar = p / q^eta

 
gen theta0 = 1/n	// cournot
gen theta1 = 1		// cartel
gen theta2 = 1/n + (1-1/n)*0.4	// pre-cartel collusion
gen theta3 = 1/n + (1-1/n)*kappa_dq4	// actual conduct
gen theta4 = 1/(n+1) + (1-1/(n+1))*kappa_dq4	// conduct with one more firm

sum theta*


* Calibrate uwm and uwk (material and capital prices per unit)

preserve 
keep pbar eta theta0 uwbar psi   bl bm bk q 
keep if pbar>0 & uwbar>0 
drop if q ==. | psi ==. | uwbar ==.  
gen N = _N 
export delimited using  uwm.txt,   replace   delimiter(tab) novar
restore 

keep if pbar>0 & uwbar>0 
drop if q ==. | psi ==. | uwbar ==.  

* don't know uwm and uwk. so need to calibrate material and capital prices. run matlab code for this.

cd /Applications/MATLAB_R2022b.app/bin/		
 
 ! ./matlab -r "run $station/uwm_master.m; "	// this matlab code calibrates uwm and uwk such that predicted and observed output is matched
    
cd "$station"


preserve 
import delimited using  uwm_hat.txt, clear 
rename (v1) (uwm)
global uwm = uwm
restore

gen uwm = $uwm 	
gen uwk = uwm
  
 
forvalues m = 0/4 {
gen mu`m' = 1/(1+eta*theta`m')
gen q`m' = (((bl*pbar*(1+eta*theta`m'))/(uwbar*(1+psi*theta`m')))^(bl/(1+psi)) *((bm*pbar*(1+eta*theta`m'))/(uwm))^(bm) *((bk*pbar*(1+eta*theta`m'))/(uwk))^(bk))^(1/(1-(1+eta)*bl/(1+psi) - (1+eta)*bm - (1+eta)*bk )) 
gen emp`m' = (bl*q`m'^(1+eta)*pbar*(1+eta*theta`m')/(uwbar*(1+psi*theta`m')))^(1/(1+psi))
gen uw`m' =  uwbar * emp`m'^psi 
gen p`m' = pbar *q`m'^eta
gen wm`m' = (bm* p`m'*q`m') / (mu`m')
gen wk`m' = (bk* p`m'*q`m') / (mu`m')
gen m`m' = wm`m'/uwm
gen k`m' = wk`m'/uwk
} 
sum pbar
sum q0  q


* estimate fixed costs 
 
forvalues m = 0/4 {
gen rev`m' = q`m'*p`m' 
gen cost`m'= emp`m'*uw`m' + wm`m' 
gen vp`m' = q`m'*p`m' - emp`m'*uw`m' - wm`m' 
}

 
sum vp*
 
gen fc_upper = vp3/n  
gen fc_below = vp4/n

 

gen fcm_upper = fc_upper/1000000
gen fcm_below = fc_below/1000000

sum fcm_upper fcm_below, d

egen medfcm_upper = median (fcm_upper)
egen medfcm_below = median (fcm_below)

estpost su medfcm_upper medfcm_below  
est store fc_en

gen fc_mid = (fc_upper + fc_below)/2
sum fc_mid, d 

sum uw0 uw1 uw2  uw 
sum emp0 emp1 emp2 emp
sum q0 q1 q2 q
sum p0 p1 p2 p


gen lfc_mid = log(fc_mid)
gen lwinv = log(winv)

sum fc_mid winv 
reg lfc_mid lwinv
corr lfc_mid lwinv // 

label var lfc_mid "log(Estimated fixed cost)" 
label var lwinv "log(Capital investment)" 

gen x45 = lfc_mid 
label var x45 "45 degree line" 

twoway scatter lfc_mid lwinv||line lfc_mid x45, graphregion(color(white)) mfcolor(white)  xtitle("log(Capital investment)") legend(off) xscale(titlegap(*5)) ylabel(,angle(horizontal)) 
graph export ./output/fig/fig_fc_winv_end.pdf, replace  
 

forvalues n = 0/4 {
gen dx`n' = fc_mid*n > vp`n'
replace dx`n' = . if fc_mid==.|vp`n'==.
}
sum dx* 					 


 
* save exit probabilities to table
 
foreach var of varlist dx* {
egen av`var'=mean(`var')	
}

 
gen ratx01 = (avdx0/avdx1  - 1)
 
gen ratx21 = (avdx2/avdx1 - 1)

sum ratx*
 
foreach var in "ratx"  {
gen `var' = `var'01 	
}

estpost su ratx 

est store ratx01_en 

foreach var in "ratx"  {
replace `var' = `var'21
}

estpost su ratx 
est store ratx21_en 

label var ratx "Relative exit change"

label var ratx "Relative exit change"
  
label var ratx "Relative exit change"

label var medfcm_upper "Upper bound"
label var medfcm_below "Lower bound"
  
esttab fc_ex fc_en  using ./output/tab/table_exit.tex, replace  ///
mtitle("Exg. price" "End. price"  ) nolines prehead(\textit{Panel A: Fixed costs }    &  \multicolumn{2}{c}{Average fixed cost (million BEF)}      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&   ///
\\     \hline    )  posthead( \hline  &&&\\   ) 

esttab ratx01_ex ratx21_ex  using ./output/tab/table_exit.tex, append  ///
mtitle("Cournot" "Pre-1898 conduct"  ) nolines prehead(\textit{Panel B: Exit change - exogenous price }    &  \multicolumn{2}{c}{Change from cartel to:}      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&   ///
\\     \hline    )  posthead( \hline  &&&\\   ) 

  
esttab ratx01_en ratx21_en  using ./output/tab/table_exit.tex, append  ///
mtitle("Cournot" "Pre-1898 conduct"  ) nolines prehead(\textit{Panel C: Exit change - endogenous price }    &  \multicolumn{2}{c}{Change from cartel to:}      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&   ///
\\     \hline    )  posthead( \hline  &&&\\   ) 


sum vp* 

gen ratw10 = uw1/uw0 -1	// collusion vs. cournot
gen ratl10 = emp1/emp0 -1	//  

gen ratw12 = uw1/uw2 -1	// collusion vs. pre-cartel
gen ratl12 = emp1/emp2 -1	//  

gen ratp10 = p1/p0 -1	// collusion vs. pre-cartel
gen ratp12 = p1/p2 -1	//  
 
gen ratq10 = q1/q0 -1	// collusion vs. pre-cartel
gen ratq12 = q1/q2 -1	//  

gen ratm10 = wm1/wm0 -1	// collusion vs. pre-cartel
gen ratm12 = wm1/wm2 -1	//  


sum rat*

foreach var in "ratw" "ratl" "ratp" "ratq" "ratm" {
gen `var' = `var'10	
}

estpost su ratw ratl ratp ratq  
est store cf10_en 

foreach var in "ratw" "ratl" "ratp" "ratq" "ratm" {
replace `var' = `var'12	
}

estpost su ratw ratl ratp ratq  
est store cf12_en

label var ratw "Relative wage change"
label var ratl "Relative employment change"
label var ratq "Relative output change"
label var ratp "Relative price change"

 
esttab cf10_en cf12_en  using ./output/tab/table_cf.tex, append ///
mtitle("Cournot" "Pre-1898 conduct"  ) nolines prehead(\textit{Panel B: Endogenous price }    &  \multicolumn{2}{c}{Comparison of cartel to: }      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&   ///
\\     \hline    )  posthead( \hline  &&&\\   ) 


collapse(mean) ratw10 ratw12 ratl10 ratl12 ratq10 ratq12 ratp10 ratp12 ratm* fc*, by(yr dq4id)
save ./data/temp/ratw_end, replace 

* Endogenous exit  (endogenous P case)
 
use ./data/temp/data_cf_temp, clear
 
merge m:1 dq4id yr using ./data/temp/ratw_end, nogen

gen vp_obs = rev - wl - wm 

drop if vp_obs==.|ratw10==.|ratl10==.|fc_mid==.

gen dx_obs = fc_mid > vp_obs 


gen vp_cf1 = rev*(1-ratq10)*(1-ratp10) - wl*(1-ratw10)*(1-ratl10) - wm*(1-ratm10) 
gen dx_cf1 = fc_mid > vp_cf1 
 
gen vp_cf2 = rev*(1-ratq10)*(1-ratp10) - wl*(1-ratw12)*(1-ratl12) - wm*(1-ratm10) 
gen dx_cf2 = fc_mid > vp_cf2 
 
foreach var of varlist dx_cf1 dx_obs dx_cf2 {
bys yr: egen av`var' = mean(`var')
}

sum dx*

gen ratx_cf1 = (avdx_cf1/avdx_obs-1) 
gen ratx_cf2 = (avdx_cf2/avdx_obs-1) 

sum ratx*
 
 use ./data/temp/data_cf_temp, clear

* Compare growth rates of markdown estimates 
sum psic_dq4
gen lpsic_dq4 = log(psic_dq4)

preserve
collapse medmd_bb1 medpsic_dq4, by(yr)

foreach var of varlist medmd_bb1 medpsic_dq4 {
gen l`var' = log(`var')
}

reg lmedmd_bb1 yr if yr>=1880
reg lmedpsic_dq4 yr if yr>=1880
 
restore 
   
   
   
 ** (4) Reallocation towards high-markdown firms

gen realmd_bb1 = agmd_bb1-avmd_bb1	// reallocation towards high-markdown firms
gen realmu_bb1 = agmu_bb1-avmu_bb1	// reallocation towards high-markdown firms

** Plot all series 
 
 local m = 1
 sort yr
replace agmd_bb`m' = . if agmd_bb`m'==0
 twoway (connect avmd_bb`m' yr, yaxis(1) ytitle("Markdown", axis(1) /*size(large)*/) lwidth(medium) msize(medium) mcolor(white) mlcolor(navy) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) lcolor(navy)) ///
(connect agmd_bb`m' yr, lwidth(medium) msize(medsmall) msymbol(diamond) mcolor(white) mlcolor(maroon) lcolor(maroon ) lpattern(solid)) ///
, graphregion(color(white))  xscale(titlegap(*5)) xtitle("Year", /*size(large)*/) legend(order(1 "Unweighted average" 2 "Weighted average") /*size(large)*/ col(2)) xlabel(, /*labsize(large)*/) 
graph export ./output/fig/fig_agmd`m'.pdf, replace  
  
sort yr

forvalues m = 1/1 { 
 twoway (connect agmd_bb`m' yr, yaxis(1) ytitle("Markdown", axis(1) /*size(large)*/) lwidth(medium) msize(medium) mcolor(white) mlcolor(navy) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) lcolor(navy)) ///
(connect medmd_bb`m' yr, lwidth(medium) msize(medsmall) msymbol(diamond) mcolor(white) mlcolor(maroon) lcolor(maroon ) lpattern(solid)) ///
, graphregion(color(white)) xline(1897, lpattern(longdash)) xscale(titlegap(*5)) xtitle("Year", /*size(large)*/) legend(order(1 "Weighted average" 2 "Median") /*size(large)*/ col(2)) xlabel(, /*labsize(large)*/) 
graph export ./output/fig/fig_md`m'.pdf, replace
  

 twoway (connect agmu_bb`m' yr, yaxis(1) ytitle("Markup", axis(1) /*size(large)*/) lwidth(medium) msize(medium) mcolor(white) mlcolor(navy) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) lcolor(navy)) ///
(connect medmu_bb`m' yr, lwidth(medium) msize(medsmall) msymbol(diamond) mcolor(white) mlcolor(maroon) lcolor(maroon ) lpattern(solid)) ///
, graphregion(color(white))  xscale(titlegap(*5)) xtitle("Year", /*size(large)*/) legend(order(1 "Weighted average" 2 "Median") /*size(large)*/ col(2)) xlabel(, /*labsize(large)*/) 
*graph export ./oligopsony_ir/fig/fig_mu`m'.pdf, replace

  
  }
 

  	 
save ./data/temp/tempap, replace 
  
 *******************************************************************************
 *******************************************************************************
 *******************************************************************************
 *******************************************************************************
 
*** Bootstrapped standard errors
  do ./dr_main_bs.do 
  
  save ./data/temp/bsestimates, replace

 

  
 
*** Figure with bootstrapped standard errors

use ./data/temp/tempap, clear

gen sd_med = $sem_kappa_dq4
gen sd_av = $se_kappa_dq4

gen sd_rts_med = $sem_kappa_rts_dq4
gen sd_rts_av = $se_kappa_rts_dq4


sum sd*

foreach mom in "av" "med" {
gen `mom'kappa_ci5 = `mom'kappa_dq4 - sd_`mom'*1.96
gen `mom'kappa_ci95 = `mom'kappa_dq4 + sd_`mom'*1.96
gen `mom'kappa_ci10 = `mom'kappa_dq4 - sd_`mom'*1.65
gen `mom'kappa_ci90 = `mom'kappa_dq4 + sd_`mom'*1.65
}
foreach mom in "av" "med" {
gen `mom'kappa_rts_ci5 = `mom'kappa_rts_dq4 - sd_rts_`mom'*1.96
gen `mom'kappa_rts_ci95 = `mom'kappa_rts_dq4 + sd_rts_`mom'*1.96
gen `mom'kappa_rts_ci10 = `mom'kappa_rts_dq4 - sd_rts_`mom'*1.65
gen `mom'kappa_rts_ci90 = `mom'kappa_rts_dq4 + sd_rts_`mom'*1.65
}

 


 twoway line medkappa_ci10 medkappa_ci90 medkappa_rts_ci10 medkappa_rts_ci90   yr , lpattern(shortdash shortdash ".-." ".-." ) lcolor(red red blue blue) || connect medkappa_dq4 medkappa_rts_dq4 yr, lpattern(solid) msize(medium) msymbol(circle square) lcolor(red blue) mfcolor(white white) mcolor(red blue) /*xline(1901, lwidth(medthick) lpattern(longdash) lcolor(orange_red))*/ graphregion(color(white)) xline(1897, lpattern(longdash)) legend(order(5 "Collusion index (RTS free)" 6 "Collusion index (RTS = 1.05)" 1 "10%-90% C.I. (RTS free)" 3 "10%-90% C.I. (RTS = 1.05)" )) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) ytitle("Median collusion index") xtitle("Year") xscale(titlegap(*5))
graph export ./output/fig/fig_medcol_se10.pdf, replace


 twoway line medkappa_ci5 medkappa_ci95 medkappa_rts_ci5 medkappa_rts_ci95   yr , lpattern(shortdash shortdash ".-." ".-." ) lcolor(red red blue blue) || connect medkappa_dq4 medkappa_rts_dq4 yr, lpattern(solid) msize(medium) msymbol(circle square) lcolor(red blue) mfcolor(white white) mcolor(red blue) /*xline(1901, lwidth(medthick) lpattern(longdash) lcolor(orange_red))*/ graphregion(color(white)) xline(1897, lpattern(longdash)) legend(order(5 "Collusion index (RTS free)" 6 "Collusion index (RTS = 1.05)" 1 "5%-95% C.I. (RTS free)" 3 "5%-95% C.I. (RTS = 1.05)" )) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) ytitle("Median collusion index") xtitle("Year") xscale(titlegap(*5))
graph export ./output/fig/fig_medcol_se.pdf, replace


 twoway line avkappa_ci5 avkappa_ci95 avkappa_rts_ci5 avkappa_rts_ci95   yr , lpattern(shortdash shortdash ".-." ".-." ) lcolor(red red blue blue) || connect avkappa_dq4 avkappa_rts_dq4 yr, lpattern(solid) msize(medium) msymbol(circle square) lcolor(red blue) mfcolor(white white) mcolor(red blue) /*xline(1901, lwidth(medthick) lpattern(longdash) lcolor(orange_red))*/ graphregion(color(white)) xline(1897, lpattern(longdash)) legend(order(5 "Collusion index (RTS free)" 6 "Collusion index (RTS = 1.05)" 1 "5%-95% C.I. (RTS free)" 3 "5%-95% C.I. (RTS = 1.05)" )) ylabel(,angle(horizontal) axis(1) /*labsize(large)*/) ytitle("Average collusion index") xtitle("Year") xscale(titlegap(*5))
graph export ./output/fig/fig_avcol_se.pdf, replace


  
 
*** Write tables

use ./data/temp/bsestimates, clear

sum _th*  
   
** Table 1: production function

foreach coef in "bl" "bk"   "blk"   "bmk"  "thl" "thm" "bm" "bd" "scale" "rho"{
gen `coef' = .
}
 
label var rho "Serial correlation TFP \hspace{2em}\hspace*{\fill}    $\rho$"
label var bl "log(Labor) \hspace{2em}\hspace*{\fill}    $\beta^l$"
label var bk "log(Capital) \hspace{2em}\hspace*{\fill}   $\beta^k$"
label var bm "log(Materials) \hspace{2em}\hspace*{\fill}    $\beta^m$"
label var blk "log(Labor)*log(Capital) \hspace{2em}\hspace*{\fill}   $\beta^{lk}$"
label var bmk "log(Materials)*log(Capital) \hspace{2em}\hspace*{\fill}   $\beta^{mk}$"
label var thl "Output elasticity of labor  \hspace{2em}\hspace*{\fill}   $\theta^{l}$  "  
label var thm "Output elasticity of materials  \hspace{2em}\hspace*{\fill}   $\theta^{m}$  "  
label var scale "Returns to scale  \hspace{2em}\hspace*{\fill}   $\nu$  "  

* Cobb-Douglas

esttab pf_ols1 pf_ols1_se pf_bb1 pf_bb1_se  pf_rts1 pf_rts1_se   using ./output/tab/table1.tex, replace ///
mtitle("Est." "S.E." "Est." "S.E." "Est." "S.E.") prehead(   \textit{Panel A: Production function  }  & \multicolumn{2}{c}{log(Output)} & \multicolumn{2}{c}{log(Output)} & \multicolumn{2}{c}{log(Output)} \\    )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs ///
prefoot( &&&&\\ Method &	\multicolumn{2}{c}{OLS} &	\multicolumn{2}{c}{GMM} &	\multicolumn{2}{c}{GMM}\\ RTS &	\multicolumn{2}{c}{Free} &	\multicolumn{2}{c}{Free} &	\multicolumn{2}{c}{Fixed at 1.05}  ///
\\  R-squared  &	\multicolumn{2}{c}{`r2_ols1'} &	\multicolumn{2}{c}{`r2_gmm1'} &	\multicolumn{2}{c}{`r2_rts1'} \\  Hansen J-test &	\multicolumn{2}{c}{ } &	\multicolumn{2}{c}{`j_gmm1'} &	\multicolumn{2}{c}{`j_rts1'} \\  Hansen J-test p-value &	\multicolumn{2}{c}{ } &	\multicolumn{2}{c}{`jp_gmm1'}&	\multicolumn{2}{c}{`jp_rts1'} \\ ///
 No. firms &	\multicolumn{2}{c}{`Nf_ols1'} &	\multicolumn{2}{c}{`Nf_gmm1'}&	\multicolumn{2}{c}{`Nf_rts1'} ///
\\  Observations &	\multicolumn{2}{c}{`N_ols1'} &	\multicolumn{2}{c}{`N_gmm1'}&	\multicolumn{2}{c}{`N_rts1'}\\ &&&&\\  \hline )  posthead( \hline  &&&\\   ) 
 
  
  
 
** Table 2:  markdowns and markups

*2a level

local m =1
foreach coef in avmu avmd agmu medmu agmd medmd { 
gen `coef' = .
}
 
label var avmu  "Average markup"
label var medmu  "Median markup"
label var avmd  "Average markdown"
label var medmd  "Median markdown"
 
esttab md1_ols md1_ols_se md1_bb md1_bb_se md1_rts md1_rts_se    using ./output/tab/table1.tex, append ///
 mtitle("Est." "S.E." "Est." "S.E." "Est." "S.E.") prehead( \textit{Panel B: Markdowns/markups }    &  \multicolumn{2}{c}{}   &  \multicolumn{2}{c}{ }  &  \multicolumn{2}{c}{ }     \\     )    ///
cells(mean(fmt(3))   ) label booktabs nonum collabels(none) gaps f   noobs ///
prefoot(&&&&\\ Method &	\multicolumn{2}{c}{OLS} &	\multicolumn{2}{c}{GMM} &	\multicolumn{2}{c}{GMM}\\ RTS &	\multicolumn{2}{c}{Free} &	\multicolumn{2}{c}{Free} &	\multicolumn{2}{c}{Fixed at 1.05}  \\  &&&&\\  \hline   ) /// 
posthead( \hline   &&\\   ) 
 
 * labor supply table
gen bw = .
 label var bw "log(Employment)"
esttab labsup_ols labsup_ols_se labsup_iv labsup_iv_se  using ./output/tab/table1.tex, append  ///
mtitle("Est." "S.E." "Est." "S.E.") nolines prehead(\textit{Panel C: Labor supply }    &  \multicolumn{2}{c}{ log(Wage)}   &  \multicolumn{2}{c}{ log(Wage)}      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&\\ Method &	\multicolumn{2}{c}{OLS} &	\multicolumn{2}{c}{IV}  ///
\\  First-stage F-statistic && &	\multicolumn{2}{c}{`F_ls_dq4'}\\Hansen J-test&	\multicolumn{2}{c}{ } &	\multicolumn{2}{c}{`j_iv_dq4'} \\Hansen J-test p-value &	\multicolumn{2}{c}{ } &	\multicolumn{2}{c}{`jp_iv_dq4'}  ///
\\   Observations &	\multicolumn{2}{c}{`N_ls_dq4'} &	\multicolumn{2}{c}{`N_ls_dq4'}  ///
\\ Firm-level elasticity &	\multicolumn{2}{c}{`supelast_ols_f_short'} &	\multicolumn{2}{c}{`supelast_iv_f_short'} \\ \hline    )  posthead( \hline  &&&\\   ) 
 drop bw

 
 
  

** Table 3: collusion

* cross-section
gen th_car = .
gen th_union = .

label var th_car "1(Cartel)"
label var th_union "1(Employers' Association)"

* time series

 label var th_lyr "Year"
 label var th_yr1 "1(1855\$<\$Year\$<\$1865) "
 label var th_yr2 "1(1865\$<\$Year\$<\$1875)  "
 label var th_yr3 "1(1875\$<\$Year\$<\$1885) "
 label var th_yr4 "1(1885\$<\$Year\$<\$1895) "
 label var th_yr5 "1(1895\$<\$Year\$<\$1905) "
 label var th_yr6 "1(1905\$<\$Year\$<\$1915) "

esttab col1 col1_se  time2 time2_se     using ./output/tab/table3.tex, replace ///
mtitle("Est." "S.E." "Est." "S.E.") prehead( \textit{Panel A: Markdown correlations}     &  \multicolumn{2}{c}{log(Markdown)}  &  \multicolumn{2}{c}{log(Markdown)}  \\       )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs ///
prefoot( &&&&\\    Year FE &	\multicolumn{2}{c}{Yes} &	\multicolumn{2}{c}{No}\\  R-squared  &	\multicolumn{2}{c}{`r2_col1'} &	\multicolumn{2}{c}{`r2_time2'} \\  Observations &	\multicolumn{2}{c}{`N_col1'} &	\multicolumn{2}{c}{`N_time2'} \\ &&&&\\  )  posthead( \hline  &&&\\   ) 

esttab col3 col3_se   col4 col4_se     using ./output/tab/table3.tex, append ///
mtitle("Est." "S.E." "Est." "S.E.") prehead(\hline \textit{\textit{Panel B: Employers' assoc.: pre- vs. post-cartel}}   &  \multicolumn{2}{c}{log(Markdown)}  &  \multicolumn{2}{c}{log(Markdown)}  \\       )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs ///
prefoot( &&&&\\   Time period &	\multicolumn{2}{c}{1845-1897} &	\multicolumn{2}{c}{1898-1913}\\     R-squared  &	\multicolumn{2}{c}{`r2_col3'} &	\multicolumn{2}{c}{`r2_col4'} \\  Observations &	\multicolumn{2}{c}{`N_col3'} &	\multicolumn{2}{c}{`N_col4'} \\ &&&&\\ \hline     )  posthead( \hline  &&&\\   ) 

 
 
* Markdown-size correlation 

gen mdcorr = .
label var mdcorr "log(Labor market share)"
 esttab mdcorr_1a mdcorr_1a_se mdcorr_1b mdcorr_1b_se mdcorr_1c mdcorr_1c_se using ./output/tab/table_mdcorr.tex, replace  ///
mtitle("Est." "S.E." "Est." "S.E." "Est." "S.E." ) nolines prehead(\textit{Panel A: All firms }    &  \multicolumn{6}{c}{log(Markdown) }      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&\\ Fixed effects &	\multicolumn{2}{c}{None} &	\multicolumn{2}{c}{Market} &	\multicolumn{2}{c}{Market  $\times$ Year} ///
\\ R-squared & 	\multicolumn{2}{c}{`r2_mdcorr1a'}& 	\multicolumn{2}{c}{`r2_mdcorr1b'} & 	\multicolumn{2}{c}{`r2_mdcorr1c'}  ///
\\ Observations & 	\multicolumn{2}{c}{`n_mdcorr1a'}& 	\multicolumn{2}{c}{`n_mdcorr1b'} & 	\multicolumn{2}{c}{`n_mdcorr1c'}  ///
\\     \hline    )  posthead( \hline  &&&\\   ) 

 esttab mdcorr_2a mdcorr_2a_se mdcorr_2b mdcorr_2b_se mdcorr_2c mdcorr_2c_se using ./output/tab/table_mdcorr.tex, append  ///
mtitle("Est." "S.E." "Est." "S.E." "Est." "S.E." ) nolines prehead(\textit{Panel B: Non-cartel firms }    &  \multicolumn{6}{c}{log(Markdown) }      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&\\ Fixed effects &	\multicolumn{2}{c}{None} &	\multicolumn{2}{c}{Market} &	\multicolumn{2}{c}{Market  $\times$ Year} ///
\\ R-squared & 	\multicolumn{2}{c}{`r2_mdcorr2a'}& 	\multicolumn{2}{c}{`r2_mdcorr2b'} & 	\multicolumn{2}{c}{`r2_mdcorr2c'}  ///
\\ Observations & 	\multicolumn{2}{c}{`n_mdcorr2a'}& 	\multicolumn{2}{c}{`n_mdcorr2b'} & 	\multicolumn{2}{c}{`n_mdcorr2c'}  ///
\\     \hline    )  posthead( \hline  &&&\\   ) 

 esttab mdcorr_3a mdcorr_3a_se mdcorr_3b mdcorr_3b_se mdcorr_3c mdcorr_3c_se using ./output/tab/table_mdcorr.tex, append  ///
mtitle("Est." "S.E." "Est." "S.E." "Est." "S.E." ) nolines prehead(\textit{Panel C: Cartel firms }    &  \multicolumn{6}{c}{log(Markdown) }      \\ )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs   ///
prefoot( &&&&\\ Fixed effects &	\multicolumn{2}{c}{None} &	\multicolumn{2}{c}{Market} &	\multicolumn{2}{c}{Market $\times$ Year} ///
\\ R-squared & 	\multicolumn{2}{c}{`r2_mdcorr3a'}& 	\multicolumn{2}{c}{`r2_mdcorr3b'} & 	\multicolumn{2}{c}{`r2_mdcorr3c'}  ///
\\ Observations & 	\multicolumn{2}{c}{`n_mdcorr3a'}& 	\multicolumn{2}{c}{`n_mdcorr3b'} & 	\multicolumn{2}{c}{`n_mdcorr3c'}  ///
\\     \hline    )  posthead( \hline  &&&\\   ) 

  
** Table appendix: Markdown drivers
 
foreach coef in "th_rail" "th_tram" "th_nm1" "th_nm2" "th_nm3" "th_urban"   ///
 "th_lib" "th_soc" "th_com" "th_oth" "th_libp" "th_socp" "th_comp" "th_othp" {
gen `coef' = .
}
 
label var th_rail "1(Railroad)"
label var th_tram "1(Tramway)"
label var th_nm1 "One firm"
label var th_nm2 "Two firms"
label var th_nm3 "Three firms"
label var th_urban "1(Urban)"
 
esttab conc1 conc1_se conc2 conc2_se    using ./output/tab/table_drivers.tex, replace ///
mtitle("Est." "S.E." "Est." "S.E.") prehead( &  \multicolumn{2}{c}{log(Markdown)}  &  \multicolumn{2}{c}{log(Markdown)}  \\     )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs ///
prefoot( &&&&\\  Mine FE &	\multicolumn{2}{c}{No} &	\multicolumn{2}{c}{Yes}\\ Year FE &	\multicolumn{2}{c}{Yes} &	\multicolumn{2}{c}{Yes} \\  R-squared  &	\multicolumn{2}{c}{`r2_conc1'} &	\multicolumn{2}{c}{`r2_conc2'} \\  Observations &	\multicolumn{2}{c}{`N_conc1'} &	\multicolumn{2}{c}{`N_conc2'} \\ &&&&\\   \hline   )  posthead( \hline  &&&\\   ) 

 ** Table appendix: Wage decomposition (not used)

foreach var in "w" "m" {
gen r2_`var' = .
}

label var r2_w "log(Wage)"
label var r2_m "log(Wage markdown)"

esttab r2_1 r2_2 r2_3     using ./output/tab/table_r2.tex, replace ///
mtitle(" R$^2$" "R$^2$ " "R$^2$ "  ) prehead(  \vspace{-1em}  \\    )    ///
cells(mean(fmt(3))  ) label booktabs nonum collabels(none) gaps f   noobs ///
prefoot( && \\ Year F.E. & X & X & X      \\  Municipality F.E. &     & X & X  \\ Firm F.E.   &  & & X \\   \hline)  posthead( \hline  && \\   ) 

 
 *******************************************************************************
 *******************************************************************************
 *******************************************************************************
 *******************************************************************************
 	

	
 // Appendix: Cost shares
use ./data/temp/data_adm_clean.dta, clear

collapse (sum) wl wm winv tc, by(yr)
egen check = rowtotal(wl wm  winv) 

order wl wm winv yr

gen sum2 = wl + winv
gen sum3 = wl + wm + winv

global varlist "wl sum2 sum3"
foreach k in $varlist {
replace `k' = `k' / 1000
}

gen zero = 0

twoway (area sum3 sum2 wl yr, color(dknavy navy*0.5 navy*0.1)) ///
	, graphregion(color(white)) legend(order(3 "Labor" 2 "Intermediate input" 1 "Investment") col(1)) xtitle("Year") ytitle("(in 1000s 1910BEF)")  ///
	xline(1868 1900, lwidth(medthick) lpattern(longdash) lcolor(orange_red)) ylabel(,angle(horizontal)) xscale(titlegap(*5))
graph export ./output/fig/fig_cost.pdf, replace
